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The thermodynamics of massless ideal gas of overlap quarks has been investigated both analyt- 
ically and numerically for both zero and nonzero baryon chemical potential. Any ^-divergence 
is shown analytically to be absent for a class of actions with nonzero chemical potential. All such 
actions are shown to violate chiral invariance. While the parameter M can be shown to be irrelevant 
in the continuum limit, as expected, it is shown numerically that the continuum limit can be reached 
with relatively coarser lattices for 1.5 < M < 1.6. Numerical limitations of the existing method of 
introduction of chemical potential are demonstrated. Finally we also show that the energy density 
for the massive overlap fermions has the correct continuum limit. 
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I. INTRODUCTION 

Lattice QCD has so far provided the most reliable theo- 
retical predictions for the thermodynamics of quarks and 
gluons important for the ongoing experiments at the Rel- 
ativistic Heavy Ion Collider (RHIC), and may continue 
to do so for those at the upcoming Large Hadron Col- 
lider(LHC). While the equation of state seems[l| to ex- 
hibit robust features as one changes the number of light 
quarks, Nf, the order of the phase transition and the 
transition temperature T c seems [2[ to depend on it cru- 
cially. Indeed, the location, and even the existence of the 
critical point in the /is — T phase diagram is expected Q 
to depend on Nf, as a result of this dependence of the 
order of the transitions on Nf . Since the transition seems 
to be associated with the restoration of the spontaneously 
broken chiral symmetry at high temperatures, it is very 
important to study it using fermions having exact chi- 
ral symmetries on the lattice. The popular choices of 
the fermions employed in simulations so far have either 
no chiral symmetry (Wilson fermions) on the lattice or 
only partial chiral symmetry (staggered fermions). For 
the latter even Nf is not well defined on the lattice and 
is typically taken to be the anticipated continuum value. 
Of course, these issues are expected to become irrelevant 
in the continuum limit of vanishing lattice spacing, i.e, 
in the limit when the number of sites in the temporal 
direction becomes very large : Nt — * oo. But they are 
likely to affect the current bunch of results obtained on 
lattices up to Nt = 8. 

In view of the experimental relevance of these issues, 
it would clearly be ideal to employ fermions with exact 
chiral symmetry on lattice for investigations of the QCD 
thermodynamics. As is well-known by now, the overlap 
fermions 4] have such good chiral properties even on the 
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lattice. The corresponding fermion operator respects chi- 
ral symmetry at the expense of being highly non-local, 
making the corresponding computations rather expen- 
sive. Advances in both algorithms and the computer 
hardware may have brought such investigations closer to 
reality today. In this paper we investigate the thermo- 
dynamics of the free overlap fermions with an aim to 
examine its continuum limit both analytically and nu- 
merically. For the above mentioned practical reasons, we 
investigate numerically whether the irrelevant parameter 
M (see below for explicit definition) can be tuned opti- 
mally to recover the continuum results on the smallest 
possible lattice size. These predictions can be used in 
full QCD simulations with such fermions to do the fi- 
nite temperature calculations faster. A way to introduc- 
tion of chemical potential in the overlap formalism was 
proposed[{| and has been recently studied^] numerically 
with M = 1. It was shown that the known canonical /j? 
divergence at zero temperature did not appear. Our an- 
alytical work shows the absence of the divergence for all 
allowed M, < M < 2, for the same class [Il| of actions 
as the staggered fermions. 

The plan of our paper is as follows. Section II deals 
with zero chemical potential case. The analytical (Nt — > 
oo) results for the energy density (equivalently pressure) 
at nonzero temperatures are derived and the numerical 
result for finite Nt are presented. Corresponding results 
for nonzero chemical potential are given in the next sec- 
tion, where the quark number susceptibility are also pre- 
sented in addition. Section IV covers the case of massive 
quarks. A summary is provided in the final section V. 



II. ZERO CHEMICAL POTENTIAL 

The overlap Dirac operator [J] has the following form 
for massless fermions on asymmetric lattice with spacing 
a and 04 in the spatial and temporal directions: 

D ov = 1 + lbsgn(i b D w ) , (1) 
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where sgn denotes the sign function and 
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is the standard Wilson-Dirac operator on the lattice but 
with a negative mass term M € (0,2). The overlap op- 
erator satisfies the Ginsparg- Wilson relation [12] and has 
exact chiral symmetry on the lattice. The corresponding 
infinitesimal chiral transformations [l3j for Nf — 1 are 



aip(l - ~D ov )j 5 



Sip = 075(1 - -D ov )ip and Sip 

(3) 

ip and ip are the usual four component fermion and an- 
tifcrmion fields. They acquire a flavour index for higher 
Nf, with corresponding modification in the transforma- 
tions above similar to that in the continuum. The ex- 
pression for energy density and pressure can be obtained 
from the partition function, Z = det D ov , obtained by 
integrating the quark-antiquark fields : 



From the (anti)periodic fermion boundary conditions in 
the (time) space directions, the discrete p^s appearing 
in the equations above are seen to have the following 
allowed values : 

2ni7r 



apj = —j^—,rij — 0, .., (N — l),j = 1, 2, 3 and 
a Pi = {2n + 1)n ,n = 0,.. 1 (N T -l) (8) 



N T 

Note that the variables hi are all real. Further, a sim- 
ple algebra shows that (h 2 + h 2 ) > for all ranges of 
interest for M, a and 04. Since the sign term in eq.© is 
thus a constant, it does not contribute to the derivative 
in eq.©: it merely provides the overall sign for the en- 
ergy density. Choosing sgn(^Jh 2 + h 2 ) = 1, the energy 
density becomes 



T 2 dhxZ(Y,T) 
~V df 
d\nZ(V,T) 



and 
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dV 
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where the spatial volume V = N 3 a 3 and the tempera- 
ture T = (Ntcl^)- 1 for an N 3 x N T lattice. We restrict 
ourselves to U = 1 here to focus on the ideal gas limit. 
Noting that the sign function for a matrix is defined in 
terms of its eigenvalues, the energy density can be writ- 
ten as 
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where the chiral nature of the eigenvalue spectrum in the 
free case was used in the last line. The eigenvalues of 
the free overlap operator in the momentum space can be 
easily worked [y, [j]] out to be 



-.//;/ j \jh 2 + hfj h 5 ± i\/K 2 
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where the variables h above are given by 
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where the summations are over all the discrete sets of 
momenta on the lattice. The derivatives in the expression 
above are seen to be 
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Similarly pressure P can be computed by taking partial 
derivative with respect to a, holding 04 constant to obtain 
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The derivatives in the expression for pressure are 
hd 
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Substituting the derivatives in eq.© and eq. ([T2|) . one 
finds the expected ideal gas equation of state e = 3P, 
valid for all values of a and 0,4. We shall therefore focus in 
the remainder only on the energy density for free overlap 
quarks on the lattice and evaluate it by setting 04 = a. 
We introduce a more compact notation for doing so: 



h 5 = 9- 
h 2 = f 

d + 2gcosu! , 



COS L) 



sin 2 uj 



h 2 + h\ = d+2gcosu , (15) 
where u> = ap^ and the functions g, f and d are given by 
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It may be noted that the <?, d and / depend only on spa- 
tial momenta pj and enable us to write down the ap^ = in- 
dependence of the energy density explicitly: 



ca 



N 3 N 

(1 — COS (jJn) 



■ COS W n ) + y/d + 2g COSW„ 

sin 2 w„ (g + cos w n ) 



d + 2g cos o>„ (d + 2<? cos w„)(/ + sin oj„ ) 

(17) 

As shown in the Appendix [5] the summation over the 
Matsubara frequencies can be carried out using the stan- 
dard contour integral techniques, resulting in the energy 
density on the lattice to be, 
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where £3, £4 terms come from from the line integrals 3 and 
4 in Figured] respectively. Their explicit Nt dependence 
indicates that they contribute to the energy density on 
the lattice. However, they do not do so in the continuum 
limit, as we shall see below. 

In order to take the continuum limit of a — > 0, we 
let N,Nt — > 00 such that T and VT 3 is kept constant. 
Each summation over momenta is replaced by an integral 
in this limit: 



1 x ^ a 

N ^ ~* 2^ 



d Pj 



(19) 



Further the integration variable u> = ap^ can be traded 
for p^, pushing the the branch points at ±7r ± icosh -1 ^ 
to infinity faster than the contours 3 and 4 are pushed. 



The line integrals, and hence the terms £3 and £4 vanish. 
Since the poles at i sinh -1 y/J scale as a in this limit, they 
continue to be enclosed in the contour at a finite P4 and 
do contribute to the energy density. This can, of course, 
be explicitly checked algebraically by taking the limit of 
the eq. (|18j) to obtain the expression for the continuum 
energy density as 
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(20) 



where E = \J p\ + p\ + p\ is the energy of the massless 
quarks. 



A. Numerical evaluation 

In this subsection we investigate the lattice energy den- 
sity of the eq. (fT?!) numerically by summing over all the 
momenta. Our aim is i) to estimate the importance of 
the terms £3 and £4 in it on lattices of practical sizes, 
and ii) to find out the role M plays on finite lattices. 
In particular, it would be good to know if there exists 
a range of the irrelevant parameter M for which the en- 
ergy density converges to that of the continuum ideal 
Fermi gas on reasonable, i.e. computationally inexpen- 
sive, lattice sizes. Since we have shown the existence of 
the continuum limit for the entire allowed range of M in 
the previous subsection, it is clear that a sufficiently fine 
lattice must eventually yield the correct result for any 
M. 

In general, the dimensionless lattice energy density will 
be of the form, 



(21) 



where the coefficients A(M) and B are the usual vacuum 
and the T A contributions, while C(M) and D(M) are fi- 
nite lattice spacing artifacts. For each value of M and 
aspect ratio, defined as ( = N/Nt, the energy density 
on the lattice was calculated as a function of the Nt- 
Clearly A is the dominant contribution and its removal 
turned out to be a tricky issue governed by the precision 
of our computations. Fitting the above form to obtain 
C or D was therefore not feasible. The zero temperature 
part of the energy density was calculated from eq. (|17p by 
taking Nt — > 00 and a large spatial extent, keeping the 
lattice spacing finite. The resulting integral over w was 
done numerically for each M to estimate the zero tem- 
perature contribution. Subtracting the zero temperature 
part from the energy density and dividing the resultant 
£ by csb gives us a ratio which we employ for further 
studies. Figure Q] displays the ratio e/esB as a function 
of Nt for M = 1.55 and various aspect ratios £. A mild 
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TABLE I: e/e ss values for diff. M for £ = 2 



iV T 


M=1.0 


1.50 


1.55 


1.60 


1.65 


4 


0.630 


1.453 


1.571 


1.697 


1.828 


6 


1.194 


1.606 


1.690 


1.792 


1.914 


8 


1.316 


1.355 


1.383 


1.431 


1.506 


10 


1.268 


1.158 


1.150 


1.156 


1.186 


12 


1.206 


1.078 


1.054 


1.036 


1.033 


14 


1.160 


1.060 


1.032 


1.004 


0.983 


16 


1.129 


1.061 


1.037 


1.008 


0.979 



dependence on £ is visible for lower values but in each 
case the curve approaches to unity by Nt = 12, sig- 
nalling the onset of continuum limit. Figure [21 exhibit 
the M-dependence of the same ratio for a fixed C = 5 for 
the range 0.4 < M < 1.65. A range of 1.5 < M < 1.6 
emerges as the favoured one because all the M-dependent 
terms are seen to be minimum there and hence the con- 
tinuum limit is reached faster. On smaller lattices with 
Nt = 4—8, the lattice results are seen to be 1.6-1.8 
times larger in this range of M. For other values of M, 
the continuum limit is seen to be approached slowly; even 
an Nt = 25 seems not enough. For larger M, we also ob- 
served oscillations as Nt changed between odd and even, 
limiting our effort to increase the M-range further. The 
values oi e/esB for Nt = 4—16 and different M are tab- 
ulated in Tables [1] [IT] for easy reference. We note that 
the C = 3, 4 results are the same as that for £ = 5 as 
seen from Figure [1] In order to estimate the size of the 
1/iVj, correction term for different values of M, the same 
ratio is plotted as a function of 1 /N T in Figure [3J From 
the plot, it is evident that the correction terms die down 
very fast for 1.50 < M < 1.60 and the continuum limit 
is reached within 2-3 % already for Nt = 12 whereas for 
M = 1 they are relevant even for Nt > 12. Of course, 
the continuum extrapolation for M — 1.0 is easier than 
for 1.50 < M < 1.60 due to the nonlinearities present for 
the latter. Note, however, that energy density for at least 
three different lattice sizes with Nt = 10, 12, 14 need to 
be computed for such an extrapolation. On the other 
hand, although the extrapolation for 1.50 < M < 1.60 
is difficult due to the complex variation seen in Figure 
[31 the deviation from the continuum value is within the 
typical accuracy range of the current lattice results, mak- 
ing it an optimal range for simulations. It should also be 
noted that the corrections for the overlap fermions for 
M ~ 1.55 for Nt < 12 are smaller than compared to the 
Wilson and the staggered caseQ as well. Ref. [1] deals 
with p/psb which we showed above to be identical to the 
(■/(■SB for the overlap ideal gas. 

Filled squares in Figure [H show the percentage average 
deviations of the ratio from unity due to lattice artifact 
terms as a function of M for large Nt values(iVr > 18). 
It shows marginal dependence on M for M < 1.2 but the 
deviation itself is about 5-6 %. For larger M, the data 
show a dip, indicating clearly that the thermodynamics 



TABLE II: e/ess values for diff. M for C = 5 



7V T 


M=1.0 


1.50 


1.55 
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0.561 
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1.041 


1.018 


1.001 


1.002 
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1.123 


1.023 


0.996 


0.969 


0.950 
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1.092 


1.025 


1.001 


0.972 


0.944 
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FIG. 1: The variation of the ratio e/esB with Nt for M = 1.55 
and C = 2 - 5. 



of free fermions favours the optimum value of M to lie 
between 1.50-1.60, with a deviation of only about 2.5 % 
or lower. 

Comparing our results with other studies of thermo- 
dynamics of free fermions done with improved actions Q 
and also with overlap fermions (M = 1) in 2-D[1(| as well 
as in 4-D0, we find that i) there are larger deviations 
in higher dimensions and ii) the oft-favoured choice of 
M = 1 favours rather poorly on finite lattices. Indeed, 
one can significantly reduce the corrections to the energy 
density of overlap fermions due to the lattice artifacts 
with a proper choice of M. 



III. NONZERO CHEMICAL POTENTIAL 

The chemical potential is usually introduced as the La- 
grange multiplier to investigate thermodynamics at con- 
stant conserved number. Constructing the relevant num- 
ber operator for the overlap Dirac fermions is not easy 
due to its nonlocalityflij and may even be not unique [ltf. 
Instead of deriving the conserved number, one may make 
an inspired guess for it such that it has the right contin- 
uum limit. One such proposal for introducing the chem- 
ical potential for the overlap operator is [f| to introduce 
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FIG. 2: The variation of the ratio t/tsB with Nt for different 
M and ( = 5. 



FIG. 4: The estimated finite lattice spacing corrections for 
the energy density and susceptibility at fi = in percentage 
as a function of M. 
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FIG. 3: The variation of the ratio t/tsB with l/Nr for differ- 
ent M > 1 and £ = 5. The plusses, crosses, stars, boxes and 
triangles denote M=1.0, 1.50, 1.55, 1.60 and 1.65 respectively. 



gered and Wilson fermions, if /i was introduced naively 
as a coefficient of the conserved number. These were 
eliminated by the the choice of the exp(±/t) factors. A 
more general way to introduce the chemical potential is, 
of course, to introduce functions K(fi) and L{fi) in place 
of the factors exp(/i) and exp(— fi) respectively such that 
K(ji) = l + fi + 0(jl 2 ) and L{(j) = l-fi + 0(ji 2 ), It was 
shown [ll| that the quadratic divergences are avoided if 
K(M • L(fL) = 1. 

Here we follow that idea and introduce chemical poten- 
tial in the overlap Dirac operator through the K and L 
factors in D\y and study where the condition to eliminate 
the quadratic divergences remains the same. Introducing 



g(A) - 

2 

K(fi)+L(fi) 



Rswh6 



— R cosh 9 , 



(23) 



it in the Dyy as one would for the usual Wilson fermions: 



D ov = 1 + -y 5 sgn(<y 5 D w (jl)) 



(22) 



where the chemical potential /i = fia^ appears only as 
multiplying factors exp(/t) and exp(— jj) to the links L/4 
and u\ respectively in eq.J^J). This, of course, renders 
~f§Dw{^) to be non-hermitian, necessitating an extension 
of the usual definition of the sign function. The natural 
choice [5| was to use the sign of the real part of the eigen- 
values of 75-Dw(A i ) in the equation above. It is important 
to note that the extended sign function it is not defined 
for purely imaginary eigenvalues. Numerical simulations 
were performed for an ideal gas of overlap fermions to 
show that the above way of introducing /x does not en- 
counter any quadratic divergences at zero temperature. 
Such divergences were known to arise [lT], EBl for stag- 



one can follow through the steps of the previous section 
to find that the free overlap operator in the momentum 
space can again be written in terms of the hi of eq.0 
but with /14 and h 5 changed to : 



h 5 = M — ^^(1 — cos(apj)) — — (1 — i?cos(a4p4 — i6)) 

(24) 



/14 = i?sin(a4P4 — i 

04 



The energy density in presence of finite chemical po- 
tential is defined as 



1 / din det D 



N 3 a 3 N T V da 4 
2 /9Zn(A+A_) 



N 3 a 3 N T V <9a 4 



(25) 



(26) 
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In the following we assume that the sign function is al- 
ways defined and is +1, as for the fj, — case. We shall 
comment on this assumption later. The energy density 
is obtained using eq. (|25[) and setting a = 0,4. 



1 — R cos(w„ — i9) 



N 3 N T ^ 

Pj ,n 



d R + 2gRcos(uj n — iQ) 



R 2 sin 2 (w n - i9)(g + i?cos(w„ - id)) 
{d R + 2gRcos(uj n - iO))(f + R 2 sin 2 (w„ - iO)) 



g + Rcos(uj n — id) + \/d R + 2gR cos(uj n — iff) 

(27) 

Note that the summand in the equation above has the 
same functional form as that in eq. (|17p . Indeed the only 
changes are: d R = f + g 2 + R 2 replaces d of eq. lfTB]) . 
lu — > lu — iO and the factor R multiplies each sine/cosine 
term. Let us therefore denote the summand in eq. (|27p as 
F(R,uj — iO). Comparing eq. ([24|) with eq.((7j), and using 
the expression for the pressure given in eq. (|12|) . one again 
finds that the equation of state e = 3P also holds in the 
presence of a chemical potential, fi. An additional new 
physical observable that can be computed is the number 
density, defined as, 



1 fdln&etD 



N 3 a 3 N T \ dfi 



(28) 



In terms of h's the previous expression can be calculated 
explicitly, 



-2i 



N 3 N 7 



^ [R sin {u> n - 



Pj,n 

gRcos(oj n ~i9)+R 2 + f 



,(d R + 2gRcos(uj n - i6))(f + R 2 sin 2 (w„ - i6)) 
[g + R cos(tj n — iff) + \J d R + 2gRcos(uj n — i 
-2i 



N 3 Ni 



J2 Fn ( r > i 



(29) 



A. Analytic Results 

In order to obtain the condition for removing the di- 
vergences, we first calculate the energy density at zero 
temperature i.e for the limit — ► 00 at finite a. The 
frequency sum j^T, in eq. (|2"T|) gets replaced by the in- 
tegral t^t dw. Subtracting the vacuum contribution 
corresponding to fi = 0, i.e., R = 1,6 = 0, the energy 
density at zero temperature is given by 

ea 4 = — / F(R,uj-i0)duj- f F{uo)du 



(30) 

Choosing the contour shown in Figure O the expression 
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FIG. 5: Contour chosen for evaluating the energy density 
and the number density for nonzero of chemical potential at 
zero temperature. The thick line indicates the Matsubara 
frequencies while the filled circles denote the poles of F(R,u). 



above can be evaluated in the complex w-plane as 



-y 



irN 3 



2™^ Res F(R,Wi) 

■T 

F(R,uj)duj 



F(R,u)dw 



TT — IV 

TT-i6 



F(R,Lu)dcu 



F(uj)duj 



(31) 



The second and fourth terms cancel since F is an even 
function which satisfies F(ir + ir/) = F(—ir + irf). Hence, 
we obtain 



TTiV 3 



Pj 

F{R,uj)duj 



V 2 

F(u)duj 



(32) 



where -ii?3 is the residue of the function F(R,u>) at the 
pole —i sinh -1 (-\/7/i?) and is given by 



Rs 



v7(g + VJTW + y/d R + 2gy/f + W) 
2^JTW{d R + 2g^JTW) 

x (g + Vf + R 2 ) ■ (33) 



Noting that the number density in eq. (|29|) has the same 
pole structure as the energy density, with only the 
residues being different in the two cases, the latter can 
also be calculated in the same way to obtain, 



-T 



ttN s 



Pj 



2vri? 4 e 



K(ji) - m 



V7 



Fjy (R, u)duj + i I FN{(jj)dui 



(34) 
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where i?4 is the residue of the function F^(uj) at the pole 
— i sinh _1 (y / 7/i?) given by, 



(.9 + y f TTW + Jd R + 2g^fTW) 



2^fTW(d R + 2g^JTW) 

x (gV&Tf + R 2 + /) . 



Ra 



It is clear from both eqs. (|32|) and (|34| that the condi- 
tion R = 1 cancels the two integrals in each of them, 
yielding the canonical forms of the Fermi surface. For 
R 7^ 1, there will in general be violations of the Fermi 
surface on the lattice. Moreover, in the continuum limit 
a — > 0, one will in general have the /Lt 2 /a 2 -divergences for 
R =/= 1 in both the energy density and the number density. 
The condition to obtain the correct continuum values of 
e = /i 4 /47T 2 and n — /i 3 /37r 2 can also be seen to be the 
expected K(fi) — L(p,) — 2/t + 0(fl 2 ). Note that the ear- 
lier work[llj on staggered fermions employed the exact 
number density on the lattice which is not the case for 
the overlap fermions here. That one obtains still identi- 
cal conditions in both the cases suggests that it is indeed 
the behaviour near the continuum limit which dictates 
these conditions. 

Another crucial difference is that the introduction of 
the functions K and L for the staggered fermions still 
leaves the action invariant under the chiral transforma- 
tions due to the locality of the action. This is true 
for the full theory, i.e., even after the link variables , 
U£ are restored. On the contrary, one can easily check 
that one breaks the chiral invariance in the case of the 
overlap fermions by these functions K 1 L 1 or exp(±/2). 
As defined in cq.([3|), the chiral transformation involves 
D ov (fi — 0), while the action for fi ^ for the overlap 
fermions has D ov (p,) of the ea.(|22p. By construction, the 
latter does satisfy the Ginsparg- Wilson relation [l2j with 
the /x-dependent overlap Dirac operator on both sides : 



{75, -Dot) (A)} = D ov (fl)jr J D ov (fl) 



(35) 



Unfortunately though, it is not sufficient to guarantee 
invariance under the chiral transformation in eq.Q, as 
it does not have any /i-dependence. Indeed, the variation 
of action under the chiral transformation of eq.([3]) is 

SS = a^ip x \j5D0v (A) + At- (A)75 

- ~A„(o)7 5 A„(A) - 5A»(A)7isA»(o)]_A(?6) 



which clearly does not vanish in spite of eg. ([55]). 

Finally, using the same techniques to evaluate the Mat- 
subara frequencies sum, the energy density at non-zero 
temperature and chemical potential can be computed an- 
alytically. Further details are given in the Appendix [B] 



The final expression is, 

N3 



N 3 

Pi 

Vf 



1 



^/TTf e^ 11111 ' 1 v7-£)-Mr + 1 

f 



V /TT/ e ( sillh_1 \/7+£)JVt + 1 

e 3;i + e 4fj] 



(37) 



The terms €3^ and £4^ are contributions of the line inte- 
grals in Figures l9l and ITOl They vanish in the continuum 
limit leaving only the contribution due to the residues of 
the poles shown in those figures. It can be shown ex- 
plicitly that the expression for the lattice energy density 
reduces to the well known [14j result in the continuum. 



(2tt) 3 



e T 



(2tt) 3 



E-1.1 

1 + e t 



(38) 



B. Numerical evaluation 

As in section II. A, our aim of presenting numerical 
results by carrying out the sums over all momenta in 
eq. (|B7|) is to find out the importance of lattice artifacts 
in form of the terms involving £3^ and £4^, resulting from 
the line integrals 3 and 4, and to look for the role of M. 
The focus here is, of course, on the chemical potential. 
We therefore consider two observables here. One is the 
change in the energy density, Ae(/i, T) = e(fi, T) — e(0, T). 
In continuum it is given by 



AeQVT) 

ji4 



47r 2 T 4 2T 2 



(39) 



The other quantity we consider is the quark number sus- 
ceptibility at (i = 0. For the free overlap fermions, it is 
given for any £1 by 



f 



X 



N 3 a 2 N T 



d 2 lndetD 



dfi 2 



(40) 



which can be worked out to be 



2 1 



N T N 3 a< 



■E 

Pi ,P4 



[h 2 hi + h^hs cos(ap4 — ifij) u 



■sH-s-h 5 y 



+ 



s 2 (s - h 5 ) 



(41) 



u and v are in the expression above are 

0/ u \(u dhi , u dh5 \ 
u = 2(s - h 5 )(h4— + h 5-QT-) 

ds dhc. 



(42) 



dh 



4 (2h 2 + h 2 + h^cos(ap4 — ifi)) 



dh 

+ 114-— cos(ap4 — iu) + i/i4/i5sin(ap4 — iu) , (43) 

OfJ, 
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and 

s 2 = h 2 + hj . (44) 
Again in the continuum, the susceptibility is known to 

be 

u 2 T 2 

IT 1 3 

Our computations for the energy density were per- 
formed keeping the ratio r = fi/T = (mNt fixed, yielding 
a constant Ae/T 4 in the continuum from eq. (|39p . Our 
choices of r were restricted by the fact that on lattices 
with odd Nt, eigenvalues of ^Dwifi) can turn purely 
imaginary for sufficiently large fi. This is related to the 
fact that (75-Dw)*75-Dw has h 2 + h\ as eigenvalues and 

Re (h 2 + h\) = g 2 + 1 + / + 2g cos u cosh/t 

Im (h 2 + h 2 ) = 2g sin us sinh/i . (46) 

has zero imaginary part at u) — tt with negative real part 
for /i > /i c . The sign function is undefined for such cases. 
Indeed, in the interacting case it may even be possible to 
get such purely imaginary argument of the sign function 
of the overlap Dirac operator for all Nt- 

From the plots of the ratio of the Ae/T 4 on the lattice 
and in the continuum (eq. (|39p h shown in Figure [5] we 
again conclude that the continuum limit is reached for 
Nt > 12 for both the cases for essentially all M, with the 
1.5 < M < 1.6 region displaying the smallest deviations 
in the region Nt > 12 as in the = case in section II. A. 
Moreover the results for Ae again appear about 1.6-1.8 
times larger on the lattices with Nt = 4 — 8 while the 
susceptibility is close to twice the continuum result. 

In the upper panel of Figure we display the ratio of 
x(p = 0, T) on the lattice with the corresponding contin- 
uum value from eq. (|45p . One sees again a similar pattern 
as for the energy density. As in eq. (p?Tj) . the susceptibility 
calculated on the lattice will also have a form, 

, n , B x C*(M) D*(M) 

^ 0) =N 2+ ^r + ^r + - (47) 

where the only difference is the absence above of a con- 
stant term like A. Keeping only the first term, one will 
again get the effective B x to become M -dependent; its 
deviation from 1/3 will be a measure of the finite lat- 
tice spacing effects. The filled circles in Figure [4] display 
these artifact effects as a function of M which were ob- 
tained by assuming a constant behaviour in the range 
18 < Nt < 32. The absence of a dominant term like A 
in the equation above allowed us to re-do the fit with the 
inclusion of the next term for each M. We found that 
the resultant B x is already M-independent and close to 
1/3 in each case. Moreover the C x changed with M sub- 
stantially and was smallest for M = 1.6. From all these 
fits, it also emerged that by Nt = 64 the contribution of 
the C x -term becomes negligible. The lower panel of Fig- 
ure [7J exhibits the results of our attempt to verify this by 
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FIG. 6: The variation of the lattice Ae/T 4 with 1/iVf. for 
ji = 0.5//Vr(upper panel) and p, — 0.8/-/Vr(lower panel). 



extending the computations to larger lattices. One does, 
indeed, find a convergence to the continuum result irre- 
spective of the value of M from lattice sizes of 320 3 x 64. 
Note that one finds very similar effects of finite lattice 
spacing for both the susceptibility and the energy den- 
sity at fi — in Figure 21 with M ~ 1.6 emerging as 
a good choice for calculations on lattices with small Nt 
due to smallest values of the correction terms. 



IV. MASSIVE OVERLAP FERMIONS 

While we restricted ourselves to the thermodynamics 
of massless overlap fermions, most of our treatment goes 
through for the massive fermions as well. In this section 
we outline this for the /x = case. For the sake of novelty, 
we use an alternative way of doing the computation. The 
overlap operator for fermions of mass m is written as 

777 777 

Dov = {1+ 2M ) + {1 ~ 2M>9nh 5 D w ) . (48) 
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FIG. 7: The variation of x(0)/xss vs l/iV|>(upper panel) and 
vs Nt (lower panel) for different M and £ = 5. 



The eigenvalues of the overlap-Dirac operator change 
from \± in eq.© to A± -*• A±(l - m/2M) + m/M. As 
a result the energy density modifies from eq. © to 



E 



where 

a = 2 i 1 - 



(49) 



4A/ 2 



and 7 = 21 



AM 2 



Substituting the values of hi, h§ and their derivatives, 
one obtains 



2a 



N 3 Nn 



E 

Pi,n 

(1 - cosw n )(/ + sin 2 ujn) 



(d + 2gcosui n )( r Y 2 (d + 2gcosuj n ) — a 2 (g + coscj„) 2 ) 
^ sin 2 oj n (g + cos uj n ) . 

(d + 2gcosa; n )(7 2 (d + 2fifcosa; n ) — a 2 (g + cosw n ) 2 ) 

x d + 2g cos w„ + a(g + cos uj n ) . (50) 



Note that setting m = 0, reduces a = 7. Substituting 
in the equation above, and using the relation d — g 2 + 
/ + 1, it becomes identical to the expression in eq. (fl~7| . 
as expected. 

One can again use the same contour method for evalu- 
ating the energy density. By comparing with eq. (|17p , the 
functions F\ and F2 can be identified as the two terms 
obtained by removing the second pair of brackets of eq. 
(ISTJf . The poles (and branch cuts) of these functions can 
be seen to be the same except that the poles defined by 
lu = ±i sinh -1 -^// are now given by 



cos(w) = y ± z 
where y and z are defined as 



V = 



9& -1) 



or 



(51) 

(52) 
(53) 



As proved in Appendix [Cl z — y > 1, making abs(cosa>) > 
1 orw purely imaginary. The pole u> — icosh _1 (j; + z) lies 
on the imaginary axis while that for lu = zcosh _1 (y — z) 
lie on parallel lines shifted by ±7r. The choice of contour 
can be made similar to that in Figure [8l allowing only 
the former to contribute to the energy density. 

Setting m = 0, it is easy to verify that this approach 
also yields precisely the result in eq. (|A7j) by selecting 
the contour as in the upper half plane of Figure [8] Its 
analog for m ^ by can be obtained by computing the 
residue at the pole defined by eq. (|5~lj) . Instead of giving 
the full expression again, we only indicate how the results 
in the continuum limit arise. The pole positions can be 
computed to be 

cosc = . 1 = (l + a2 ^ 2 +m2) ), (54) 
and denoting by m' = m(M — 2)/M 

cos. = Q2 = -(l + a2 ^ 2 +m,2) ), (55) 

where p 2 = pf +p% +p 2 - The pole at a\ has, at order a, 
the residue 



Res Fi(ai) = ± 



a\J 'p 2 + m 2 



Note that all the other poles, including the poles at 0:2, 
and the branch cuts do not contribute to the contour in- 
tegrals, as seen in Figure[U Therefore the energy density 
in the continuum comes out to be the same as in eq. (|20|) 
but with E = \J ni 2 + p\ + p\ + p\. 



V. SUMMARY 

Investigating the thermodynamics of QCD on lattice 
with fermions which possess both the chiral symmetry 
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and the flavour symmetry relevant to our world has im- 
portant consequences for both the experimental aspects 
of the heavy ion collisions and the theoretical aspects of 
the fi — T phase diagram. Staggered fermions used in the 
bulk of the work so far are not adequate to resolve some 
of these issues. Overlap fermions, while computationally 
more expensive, may prove better in such studies in near 
future. 

We have presented analytical and numerical results on 
the the thermodynamics of free overlap fermions in 4-D 
both for zero and nonzero (baryonic) chemical potential 
by varying the irrelevant parameter M. From the energy 
density computed on the lattice in these cases, we showed 
that the expected continuum limit is reached. General- 
izing the proposed action Q for nonzero /i, we demon- 
strated that the /i 2 -divergence in the continuum limit 
is avoided for a class of functions K(fi) and L(fi) with 
K (/t) • L(fi) = 1; the choice exp(±/t) for K,L also sat- 
isfies this condition is therefore shown to be free of any 
/i 2 -divergence. In all these case, however, the chiral in- 
varinace of the action is lost for nonzero //. 

While the sign function in the free overlap Dirac op- 
erator remains a constant in computations for /i = 0, we 
pointed out that it becomes undefined on lattices with 
odd number of temporal sites for /i > /i c , where the value 
of /i c depends on M. Our numerical computations were 
restricted to smaller [i- values. The numerical results were 
mildly dependent on the aspect ratio of the spatial and 
temporal direction but changed significantly as a func- 
tion of the irrelevant parameter M of the overlap Dirac 
operator. For the choice of 1.5 < M < 1.6, both the 
energy density and the quark number susceptibility com- 
puted for fi — exhibited the smallest deviations from 
the ideal gas limit, as seen in Figure [U As seen from 
Figures [U O and lattice results approximate the con- 
tinuum well for lattices with 12 or more temporal sites, 
with typically a factor sa 1.8 larger results for smaller 
lattices with 6-8 temporal sites. It would be interesting 
to check whether the optimum M-range is still the same 
in the presence of gauge fields. 
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APPENDIX A: ENERGY DENSITY FOR n = 

Before we show the details of the energy density cal- 
culation, let us list certain useful relations amongst the 
quantities g,d,b and c introduced in eq. (|16|) . (which will 
be useful for the calculations below) : 

• Since cos(apj) < 1 for any j, g < for M < 1, 



>. 



ico.sh 



ismlr 1 7/ 



X X X X 



-ismlr 1 ~fs 



-icosir 1 f 



Imw— 



FIG. 8: The choice of contours for evaluating the tj-suni in 
cq, (117[) . The dashed lines represent branch cuts. The crosses 
denote the Matsubara frequencies uj n , while the filled circles 
denote the poles of F(uj). 



• g z + I + i 

squares, 



d > , since f is a sum of 



d 2 /4 9 2 - / - 1 = (.g 2 - / - 1)72<? 2 

=► > i + /, 

cosh -1 7^ > sinh -1 -*//- This follows trivially from 
the line above. Since > (1 + /), it follows ^ > 



2« 



VT+f ( £ < -VTT 

cosh(sinh -1 y / /) 



J, for g < 0). Noting that 
-f /, one has cosh -1 ^) > 



sinh _1 v7 (cosbT 1 ^) < -sinh^V/ for g < 0) 

The last line justifies the drawing of the contour in Figure 
[8] by avoiding the poles/cuts at ±i cosh _1 (^). 

The summation over all Matsubara frequencies oj n = 
(2k+i)7t can k e ^ Qne usm g standard textbookflij 
metnod of contours in the complex u-plane. For a general 
function F(u>), which in addition may, and does, depend 
on other variables such as pj, but this dependence will 
not be shown explicitly below, the frequency sum there- 
fore is, 



F{Lo)duj 



F(m)duj ' r " 

p iujN T 4- 1 ' / piuNT J- 1 



2tt 

At 

( A1 ) 

where the integrals are evaluated on the contour lines 
running parallel to the real axis. The second integral can 
further be re-written as 



7r-He 



— 7r-He 



F{uj)duj 
e lu}N T + 1 



F^duj- I — ^ 



7T + 2E ° > 1 

(A2) 
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The summand F(uj) in eq. (fl7|) can be split in to two 
terms, 



with 

F x (w) = 



F(u) = Fi(w) + F 2 (co) , 
(1 — cosw) sin 2 uj(g + cos a;) 



d + 2gcosui (d + 2gcosw)(/ + sin 2 uj) , 
x (g + cosuj) , (A3) 



and 

F 2 (w) 



(1 — cos a;) 



sin 2 w(p + cos a;) 



d + 2g cos a; (d + 2g cos uj)(f + sin 



x yd + 2g cos w . 



(A4) 



Both the functions Fi have a finite number of poles at 
lu = ±i sinh -1 -^// an d ±m7r ± i sinh -1 \/7, where m 
is an integer. Furthermore, Fi has poles for ^ > 

at lu — ±kn ± i cosh -1 ^- while F 2 has branch points 
at the same locations. Similarly, for ^ < the poles 
(branch points) of F\ (F 2 ) are at ±icosh -1 ^. In the 
rest of the complex u> plane both the functions are ana- 
lytic. In view of these properties the contours in eq. (|Al[) 
can be deformed to the contours shown in Figure [5] 
We chose each contour such that it lies below (above) 
the cut in the upper (lower) half of the plane. As 
shown above, cosh -1 ^ > sinh -1 ^//. Defining there- 
fore 2?7 = cosh -1 7^ — sinh -1 ^// with r\ > 0, the lines 3 

and 4 are drawn through the points =f(i sinh -1 -\/7 + iff) 
respectively to avoid the cuts shown in Figure [H 

Consequently, the frequency sum in eq. (|Al[) becomes 



2tt 



F{u n ) 



-2m 



2m £ 



Res F(to) 

Im u)>0 

Res F{uj) 



Im u)<0 

F{uj)duj 



e iu>N T + i 



F(u>)du} 



e iwN T + l ^ e -iuN T + 1 



F(uj)duj 



(A5) 



Note that the line integrals along the vertical lines 
through 7r and — 7r cancel each other due to the periodic- 
ity of the function F(lo). Indeed, in general for any func- 
tion G(lo) satisfying the property, G(7r+i?;) = G(—n+ir]), 
the sum of integrals of G{lu) along opposite vertical paths 
of equal length through — 7r and 7r is identically zero. 

The residues of the function Fi^(w) at the poles lo = 
±i sinh - yfj are ±iR\ and ±iR 2 respectively where 



i?i — R 2 



v7 



(A6) 



Our choice of the contour ensures that the poles at 
±7r±i sinh -1 ^// do not contribute to the energy density. 
By taking the limit Nt ->ooon the lattice, on finds that 
the last term of eq. (|A5[) gives the quartically divergent 
vacuum contribution in the continuum limit. Defining 
the physical energy density by subtracting it off, we have 



ea 



-y 

N 3 

£3 + 64 



(A7) 



where 63, terms come from from the line integrals 3 
and 4 respectively. 



APPENDIX B: ENERGY DENSITY ATT/0 
AND n ^ 

For evaluating the energy density, we revert back to 
the choice K(p) = exp(/i) and L(p) = exp(— p) which 
has R = 1 and = p. The physical part of the energy 
density on the lattice is calculated by subtracting off the 
p = , T = contribution, 



ea 4 = 



N 3 



Pi 



N- T ^ F ^ 



* J-K 



F(uj)du 



(Bl) 

The further evaluation of the energy density can be done, 
noting that compared to the Figure [8l the Matsubara 
frequencies (2n+ 1)w/Nt in this case are displaced along 
the imaginary axis by ip in the lower half plane with 
the choice of the function l/(exp[i(u> + ip)Nr] + 1)- The 
frequency sum can be replaced by line integrals as, 



2n 



\ " 



F(w n - ip) 



—Tr—ie—ifj, 
TT — ie—ip, 

7T + /C. — / fl 



F(uj)duj 



e i(w+ip,)N T -|- \ 

F(no)duj 



p i(ui+ijl)N T 1 1 



(B2) 



The choice of the contour will be analogous to that in 
Figure [SJ but will depend on whether the pole to = 
—i sinh -1 vT is above or below the Imw = —p line. 
Therefore the frequency sum can be split into two terms, 



F(uj n - ip) 



[F< (o<j n - ip) + F> (w n - ifi)] 



(B3) 

where F > and F< are the functions with sinh 1 y/J < 
ft (sinh -1 yff > p) respectively. We have taken p < 
cosh -1 ^ because we expect that in the continuum limit 
the p to scale as the lattice spacing whereas the second 
term to tend to infinity. For the case p < sinh -1 ^//, the 
contour was chosen as shown in Figure [9] 

Using the standard trick of rewriting the second in- 
tegral in eq JB2[) as done in eq. (|A2[) the frequency sum 
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FIG. 9: The contours chosen for evaluation of the u>— sum 
for p < sinh -1 y/J. Same notation as in Figure [8] 
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FIG. 10: The contours chosen for evaluation of the u>— sum 
for jl > sinh -1 y/J. Same notation as in Figure [8] 



becomes 
2tt 



— ^2 F < ( w ™ _ *A) = ©(sinh 1 y/f - p) 



\ /nw>0 



Res F(u 



-i(u]+ip.)N T _j_ ^ 



27TZ £ 



Res 



Imuj < 
e i{u)+ip,)N T _|_ J 



gi(a;+i/i)A r T _|_ J 



■IT 

F(uj)duj 



(B4) 



The expression can be evaluated by substituting the val- 
ues of the residues of the function F(u), which are the 
sum of residues calculated in eq. (|A6j) . The integrals 3 
and 4 are along the lines Im u> = =F(sinh -1 \ff + rf). Rep- 
resenting them as £3 M ,4 M respectively, we get 



2tt — 
J^^2 F < ( w « - *A) = ©(sinh" 1 y/f - fi) 

Air Rx 



+ 



gfsinh- 1 y/J+fi)N T + 1 

4tt Ri 



gCsinh- 1 ^J-ti)N T + I 

+ / F{uj)duj + e 3M + e 4M 



(B5) 



Similarly for p > sinh -1 y/f the frequency sum is re- 
placed by integrals along the contour as shown in Figure 
[lTfl There are no poles below the Im ui = —p, line so the 
first integral in eq. (|B2[) can be replaced by a line integral 
along the line 3. Following the same steps as discussed 
for the above case, the frequency sum reduces to 



2tt 



ip) 



eip-smh- 1 ^/]) 

( 4tt Ri 

^(sinh- 1 V7-V)N T + 1 

4tt Ri 

gtsinh- 1 y/J+/i)N T _|_ 1 



F(uj)duj + e 3M + e 4 Ab6) 

-7T / 



Finally, the energy density on the lattice is obtained from 
eq. (IBl[) by substituting in the eq. (|B3[) the frequency sums 
calculated above: 



-y 

N 3 

Pj 

V7 



i 



1 



^TTJ eCsmh" 1 Vf+fi)N T + 1 

APPENDIX C 

Here we prove the claim made in eq. (|5T 
• y > since 7 > a. 



(B7) 
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• Let £ = 2p — 1, where £ > 0. z — y > 1. 

• A little algebra shows that z 2 — (y + l) 2 = £[(</ — 
l) 2 + /] + />0 which in turn implies the relation 
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